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ABSTRACT 

We present a study of cooling in radiative shocks simulated with smoothed particle hydro- 
dynamics (SPH) and adaptive mesh refinement codes. We obtain a similarity solution for a 
shock- tube problem in the presence of radiative cooling, and test how well the solution is 
reproduced in GADGET and FLASH. Shock broadening governed by the details of the nu- 
merical scheme (artificial viscosity or Riemann solvers) leads to potentially significant over- 
cooling in both codes. We interpret our findings in terms of a resolution criterion, and apply 
it to realistic simulations of cosmological accretion shocks onto galaxy haloes, cold accretion 
and thermal feedback from supernovae or active galactic nuclei. To avoid numerical overcool- 
ing of accretion shocks onto haloes that should develop a hot corona requires a particle or 
cell mass resolution of 1O 6 M , which is within reach of current state-of-the-art simulations. 
At this mass resolution, thermal feedback in the interstellar medium of a galaxy requires 
temperatures of supernova or AGN driven bubbles to be in excess of 10 7 K at densities of 
nn = 1.0 cm -3 , in order to avoid spurious suppression of the feedback by numerical over- 
cooling. 

Key words: shock waves, hydrodynamics, galaxies: formation, methods: numerical, galaxies: 
ISM, 



1 INTRODUCTION 

Radiative cooling and shocks are two important ingredients in 
galaxy formation theory ( White & Rees 1978). Whilst most codes 
used in astrophysics have facilities for handling both of these, the 
scales at which these operate and interact are challenging. We will 
start with a short tour of the processes and of the numerical codes 
we will use to simulate them. We describe a 1 -dimensional model 
problem of a radiatively cooling shock with an analytic solution 
which we model with two popular codes in astrophysical simula- 
tions, FLASH ( [Fryxell et al.|2000| an a daptive mesh refinement 
(AMR) code, and GADGET (|Springell[2005]>, a smoothed parti- 
cle hydrodynamics code(SPH; Gingold & Monag han|1977[ |Lucy| 
1977). The results of our simulations give appropriate criteria with 
which we can analyse the efficacy of our numerical schemes in 
a wide variety of astrophysical environments. We also investigate 
the mitigating factors such as the ratio of the cooling to dynamical 
times, which may enable a simulation to give approximately correct 
results when otherwise we would consider there to be insufficient 
resolution. 

There has been considerable discussion of the treatment of 
discontinuities in SPH ( |Pncel|2008l |Read et al.||2010] ), motivated 
by problems illustrated by Agertz et al. ( 2007), but the issues high- 
lighted in this paper are of a different nature. Those papers focus 
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on spurious forces introduced by the SPH scheme whilst we focus 
on evaluating the errors introduced as we approach the resolution 
limit, which are to some extent unavoidable. 



1.1 Astrophysical shocks 

The science of astrophysics is an ideal domain for the investigation 
of shock fronts on a variety of scales. Stellar winds form shocks as 
they push in to the interstellar medium. On larger scales SNe form 
very high Mach number shocks as they plough into the surrounding 
gas and form remnants. On larger scales still, galactic winds, driven 
by starbursts and active galactic nuclei (AGN), shock against the in- 
ter galactic medium (IGM). In the context of galaxy formation we 
can also consider accretion shocks, where gravitationally acceler- 
ated infalling gas shocks to form a hot corona in the dark matter 
potential well. 

Of particular interest to us in this paper are radiatively cool- 
ing shocks. To an extent, all the aforementioned shocks have radia- 
tive cooling, however the cosmological accretion shocks and SNe 
are particularly topical. In galaxy formation simulations the SNe 
at early times form remnants well below the resolution of current 
simulations and need to be modelled with subgrid physics (see Kay 
|et al.|2002| for a review of feedback methods). The cooling of the 
hot gas causes a transition from a thermally-driven to a momentum- 
driven phase, losing a significant fraction of the SNe energy. A sim- 
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ilar transition is thought to occur in the thermal to momentum tran- 
sition of winds powered by an AGN ( Booth & Schaye 2009 ). 

Cooling in accretion shocks may affect the fuelling of star for- 
mation in the host galaxy. If the gas is shocked to too high a temper- 
ature it will not cool over a Hubble time, preventing star formation 
(though non-spherical geometries may allow the gas to compress 
first and thus cool faster, see e.g. Birnboim & Dekel 2003). In a cos- 
mological simulation, however, the resolution around these cooling 
regions may be so coarse as to resolve these cooling regions with 
only a few particles. In this paper we intend to probe the effect 
of limited numerical resolution in these cases, and how these may 
affect the outcome of the simulation. 



1.2 Physical shock fronts 

Before we concentrate on the numerical aspect of cooling in 
shocks, we begin by briefly considering the processes that occur 
in a real physical shock front. A shock front is a region where 
one of the usually conserved fluid properties, entropy, is allowed 
to change. It is worth considering why such a property is otherwise 
treated as a constant, and why shocks are a special case. 

In the kinetic theory of gases, a gas is described as a large 
number of particles (e.g. atoms, molecules, ions) in constant ran- 
dom motion. The frequency of collisions defines a timescale, and 
also a typical length between collisions, the mean free path. If all 
processes acting on the gas happen on timescales much greater than 
the time between collisions, then the classical theory of adiabatics 
tells us that there will be another conserved property which, for an 
ideal gas, is pip 1 . Here, p, p and 7 are the pressure, density and 
adiabatic index, respectively. This property is a function of the en- 
tropy, and in astrophysics is often used as a proxy. 

It is worth recalling that processes which change the fluid en- 
tropy (e.g. shocks, radiative absorption, thermal conduction) will 
occur on timescales on the order of, or shorter than, the period 
between collisions (or over lengths on the order of the collision 
length). Mechanisms which heat the gas on slower time scales will 
be adiabatic processes, and will alter the thermal energy with only 
very small increases in entrop^ 

Now we come to shock fronts. A canonical example of a shock 
front would be a 1 dimensional system where the upstream fluid 
travels supersonically with respect to the down wind fluid (i.e faster 
than the thermal velocities of the particles), until it reaches the 
shock, where the majority of its mechanical energy (the bulk mo- 
tion of the particles) is converted into thermal energy. This happens 
because the pairs of particles that collide can have very different 
velocities: particles in the shock front change their energy on a 
timescale on the order of the collision time between a pair of up 
and down wind particles, which is much shorter than that between 
two down wind, or two upwind, particles. From this description we 
immediately see that physical shocks must occur over length scales 
on the order of the mean free path, which is usually much smaller 
than other physical length scales in the problem. 

The mean free path (Ax) depends upon the number density of 
particles (n) and their collisional cross section (a), as 



1 One can of course construct systems in which the time scale of inter- 
est is long enough such that viscosity, thermal diffusion, etc. dominate the 
large scales too. These problems, however, have low Reynolds and Peclet 
numbers respectively, and are the exception rather than the norm in compu- 
tational astrophysics 
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In the case of a partially or fully ionized gas, particles may interact 
on a shorter length scale (Zel'Dovich & Raizer 1967), that of the 
plasma skin depth/plasma oscillation length 
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Since the particles are not interacting via Coulomb collisions this 
is known as a 'collisionless shock'; the mechanism of interac- 
tion is the plasma oscillation (coupling together charged particles). 
Care must be taken, however, as the post-shock gas may be out of 
thermal and ionizational equilibrium for the problem in question 
(something that would not be a concern if the collision length is 
small), making these cases challenging to simulate. 

The trapping of relativistic ions between magnetic fields in 
the up and down-stream phases and subsequent acceleration is also 
believed to be the origin of the power law spectrum of high-energy 
cosmic rays, a Fermi acceleration process. 

Finally, we should complete this discussion by mentioning tur- 
bulence as a source of entropy. In general the bulk oscillations of 
fluids will occur on scales much larger than the mean free path and 
is thus unable to change the entropy. Transfer of spectral energy to 
shorter wavelengths, however, implies that eventually bulk oscilla- 
tions reach the scale of the mean free path and will be dissipated 
into thermal energy ( Kolmogo rov| 1 94 T) . 

1.3 Shocks in simulations, artificial viscosity 

Now let us consider shocks in simulations. Almost exclusively, 
the resolution of simulations will be much coarser than a phys- 
ical shock width. This is not necessarily a problem, however, as 
the bulk properties of the post- shock gas may be deduced from the 
conservation of energy and momentum, and the assumption that the 
shock process does not produce oscillations on scales much larger 
than the mean free path. 

In this paper we will contrast two schemes for numerical 
hydrodynamics that are popular in cosmology: SPH and AMR. 
Smoothed Particle Hydrodynamics (S PH; [ Gingold & M onaghan 
|T977t ; [Lucy] (1977) , see |Monaghan] (2005) and |Springel| 12005) 
for recent reviews) is a (pseudo) Lagrangian scheme in which the 
fluid is represented by a set of particles that move along with the 
flow. In this paper we will illustrate the behaviour of SPH using 
the GADGET SPH implementation of |Springel|{2005) . Adaptive 
Mesh Refinement (AMR) follows how fluid flows across a (station- 
ary) computational mesh, whose cell size may be locally 'refined' 
or 'de-refined' based on some criterion. In this paper we use the 
FLASH code, a block-structured AMR implementation by |Fryx-| 
|ell et al.|j2000l >. 

The physical process of kinetic energy dissipation by particle 
collisions is represented in the continuum approximation by a vis- 
cous term in the Navier-Stokes equations, 
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is the viscous stress tensor and 77 and ( are known as the shear and 
bulk viscosity coefficients, respectively. These coefficients can be 
measured for real fluids, however in most astrophysical flows they 
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are so small that the viscous term is insignificant outside of shocks 
(i.e. the flows have high Reynolds number). 

The variant of SPH used in this paper handles shocks with a 
prescription known as artificial viscosity (although Godunov type 
methods for SPH also exist, see Inutsuka 1994). Artificial viscosity 
was originally developed for grid codes ( |von Neumann & Richt-| 
|myer|195"0| ), and use the bulk viscosity term in Eq. ([3]), however, 
with the coefficient raised by several orders of magnitude. These 
larger values prevent the shocks generating large unphysical oscil- 
lations due to the coarseness of the sampling (see the numerical 
stability criterion of Friedrichs & Lax 1971). In SPH they also ful- 
fil a second role of preventing particle interpenetration (see |Bate| 
1 1995 1 for a thorough discussion). A number of artificial viscosity 
prescriptions are in use, the most common being that of Monaghan 
& Balsara (Balsara 1995), a Lax-Friedrichs style viscosity that is 
turned on for compressing flows. The implementation in our ver- 
sion of Gadget is based on signal velocities ( Monaghan 1997). 

In mesh codes shocks can be treated with artificial viscosity 
but more commonly a conservative Riemann solver (based upon 
Godunov's scheme, [Godunov & Ryaben ki|1964|> is used. Riemann 
solvers (see e.g the HLL solver, Harten et al. 1983) give exact so- 
lutions in Id or planar shock problems with homogeneous pre- 
and post-shock fluids, but are somewhat diffusive in other cases. 
They are still the preferred method for grid codes, however, and 
the default used in FLASH is a directionally split Riemann solver 
(Colella & Woodward 1984). Oscillations near the discontinuities 
are controlled with a monotonicity constraint. 



1.4 Radiative cooling 

Radiative cooling is an essential ingredient in galaxy formation as 
it is the process which allows the baryons in dark matter halos to 
dissipate thermal energy and thus collapse to form galaxies. Multi- 
ple cooling mechanisms are important in the astrophysical domain, 
however in this paper we will primarily be interested in collisional 
line cooling and at higher temperatures, thermal bremsstrahlung. 
The evolution of the specific thermal energy, u, due to cooling can 
be written as 



pu\ A = -A(T;Z)n 2 



(5) 



where p is the density, T the temperature, Z the metallicity and n 
the particle number density (for brevity we will subsequently refer 
to the radiative component of the specific cooling rate u\ A as it a). 
When baryon-photon interactions with the CMB and an ionizing 
background are important we have followed the prescriptions of 
Wiers ma et al.| ( |2009a] ) (see also Fig. [7] below). 

The implementation of cooling in our versions of GADGET 
and FLASH is performed by an adaptive time step integration over 
each cell/particle. The effects of cooling are included in the hydro- 
dynamic solver by operator-splitting, i.e. the separation of the two 
processes A (radiative cooling), and B (shock heating) into indi- 
vidual steps, 
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where the errors on the latter term on the order of the time step At 
depend upon the commutator [A, B]. Since the physical process of 
shock heating should occur over a much shorter time scale than that 
of radiative cooling we can justify this being zero. The numerical 
implementation of shock heating will of course take a longer time 
scale and thus would interact with the cooling if operator splitting 



was not introduced, however there is no physical motivation to pre- 
fer such a scheme in this case. 



2 RADIATIVELY COOLING SHOCKS, A MODEL 
PROBLEM 

A typical test problem in numerical hydrodynamics is that of the 
formation of a 1 -dimensional shock in a 'test tube'. In one form 
of this problem a tube is initialised with gas of constant polytropic 
index 7, the left and right halves converging with opposing veloc- 
ities vo and — vq. For a sufficiently high velocity vo a shock will 
form, creating a downstream region with higher temperature, pres- 
sure and density. This problem has a similarity solution for constant 
7. In our set-up the gas is allowed to cool radiatively, and the down- 
stream region can then cool to form a dense post- shock region. The 
initial conditions are thus 
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We note that the symmetry of this problem makes it equivalent to 
the wall shock (where there is an immovable boundary at x — 0, 
Monaghan 1997). In order to minimise the amount of modification 
in our SPH code we chose to set up the symmetric problem rather 
than implement an immovable boundary. 



2.1 Similarity solution for a radiative ID shock 

If the temperature dependence of the cooling rate is sufficiently 
simple, then this 1 dimensional shock problem has a similarity so- 
lution even in the presence of cooling. Such is the case for a cooling 
function which is a piecewise linear function of the temperature, 
such that the rate of radiative cooling, u\ a, of the specific energy u, 
is given by a 'cooling spike' : 
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where A is a positive constant and cooling is maximum at T — 
\ (Ti + To). For simplicity, in all simulations discussed below we 
initialise the temperature to To (where the cooling vanishes), so the 
initial gas is not cooling. 

The gas is chosen to be pure atomic hydrogen, i.e. 
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For comparison the reader should see the simulations of 
|Hutchings & Th omas (2000) who used a more realistic astrophys- 
ical cooling function, at the expense of not having an analytic so- 
lution. For the cooling post-shock region we find analytic solutions 
in a similarity variable A = p/ po of the form (see Appendix 
details) 
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Figure 1. Solid lines represent the analytic solution for the colliding gas problem discussed in Section [2] when cooling is included. Incoming gas from the left 
and right shocks and then compresses and cools to form a cold dense region in the center. For the example shown, the Mach number of the upstream gas is 1.5 
w.r.t. the cold, central gas and the time is 5.1 Axa/co (see Eq. fTT)). For comparison, dashed lines show the solution without cooling at the corresponding 
time. At early times (i.e. t < Ax\/co) the cooling profile is not of given form, as it has not had sufficient time to reach a stationary state (details of the 
similarity solution for cooling through a shock can be found in appendix jAl) . 
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The value of the shock velocity v s and the final density in the 
cold, post- shock region p can be found by imposing conservation 
of mass and momentum (see Appendix \A2\ . The solution is shown 
in Figure [T] 

2.2 Shock stability 



|Chevalier & Imamura| {1982) find that positive increasing linear 
cooling functions produce stable shocks. Applying this to the cool- 
ing function in Eq. fl2| ) we see that we have stable shocks pro- 



vided the post shock temperature T s < §(Ti + T ) or T s > Ti, 
which is the case for all the shocks we study later (we define the 
shock temperature T s as the temperature immediately after the 
shock, which is computed in Appendix \A2\ . If the gas is shocked 
to I (Ti + T ) < T s < Ti then the shock may be unstable as the 
cooling function has a negative slope, dr (—pu\) < 0. Intuitively 
this can be understood in terms of the length of the cooling region: 
if the length increases the shock velocity will be higher, causing 
the post shock gas to be hotter, which increases the cooling time, 
which feeds back into a longer cooling region. 



2.3 Numerical solution 

2.3. 1 Initial conditions 

The similarity solution is described with two (dimensionless) pa- 
rameters. The first is the ratio of the upper to the lower temperature 
in the cooling spike, which we will set to 20, i.e. T\ — 20T<d. This 
is motivated by the temperature dependence of the radiative cool- 
ing function of an astrophysical plasma (see also Fig. [7}, where 
individual elements contribute significantly to the cooling over ap- 
proximately a 1 dex range in temperature. 
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The second parameter is the Mach number of the shock, which 
we will quote in the rest frame of the problem (rather than the rest- 
frame of the post-shock gas, for example), 



M 



CO 



(16) 



where Co = (iPo/po) 1 ^ 2 is the upstream sound speed. Our tests 
are performed at M — 4.70 and M = 6.04. The former has been 
chosen such that the shock reaches a temperature somewhat be- 
low the maximum of the cooling function, (Ti + To)/ 2 (where the 
shock will be stable) and the latter such that the shock reaches a 
temperature somewhat above Ti (where there is no cooling). 
We plot positions in units of the cooling length, 

ksToCo 
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Similarly we express times in units of Axa/co. As observed in 
Monaghan] ( |1997| ), numerical schemes (including both SPH and 
AMR) usually produce a transient unphysical thermal bump at 
t — when there is no post- shock region. To avoid contamina- 
tion by this transient we run our simulation for a time 14.2 Ax a / Co 
and 7.1Axa/co for the M = 4.7,6.04 shocks respectively (i.e. 
we simulate for several sound crossing times of the cooling region, 
to make sure it is in a stationary state). 

For our SPH simulations we set up a long box, periodic along 
all boundaries. The particle mass is chosen to be 



mspH = Po • (0.3Ax A ) 



(18) 



(i.e a mean inter-particle spacing of 0.3 Ax a). We note that this 
setup creates a rarefaction wave that propagates inwards from the 
far edges of the computational volume (due to the discontinuity 
on this boundary), and thus we need a box long enough such that 
these rarefactions do not reach our domain of interest in the sim- 
ulation time. The particles were set up with a 'glass' distribution 
(White 1994) to minimise relaxation effects in the pre-shock fluid 
(the SPH kernel allows a cubic lattice arrangement of particles to 
slightly reduce its density, and hence release some thermal energy, 
by relaxing to a glass like state). We also raised the level of the bulk 
artificial viscosity constant, a, to 3 (from 1, see Springel 2005] for 
a complete description of the artificial viscosity prescription). We 
found this to be necessary to prevent ringing and the appearance 
of large scatter in the entropy of SPH particles in the post- shock 
region (see also Abel 2010). 

For the AMR simulations we again set up a long box with cell 
spacing 0.3Axa, with periodic boundaries in the y and z direc- 
tions and inflowing gas along the (long) x axis. No refinement was 
allowed, effectively making this a uniform Eulerian mesh. 

We considered allowing an alternative refinement criterion, 
however the standard FLASH refinement schemes will refine a 
shock to the maximum allowed level (since it contains a disconti- 
nuity), reducing it to the uniform mesh case. We refer the reader to 
the dashed lines in Figs.|2]&|4]to compare resolutions. 

We note that the use of inflowing boundary conditions in 
FLASH allowed us to avoid the rarefaction waves we created in 
SPH, and thus we could use a much shorter box (by a factor 10). 
To set the scene we begin by looking at shocks in the absence of 
cooling. 



2.3.2 Test without cooling 

The test problems in the absence of cooling are compared in the up- 
per panels of Fig. [2] and [4] (M = 4.7, 6.04 respectively). Provided 
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Figure 2. Upper panel plots the temperature in a Mach M = 4.7 shock 
without cooling. Each blue cross represents a column of FLASH cells (the 
tube is orthogonal to the mesh), each red point represents a GADGET parti- 
cle, the black line is the analytic solution. Red dashes denote the smoothing 
lengths of the GADGET particles, blue dashes the separation of FLASH 
cells (right axes). Incoming gas from the left (and right, not shown) col- 
lides to form a homogeneous hot, rarefied region in the centre. As expected, 
both codes reproduce the correct profile relatively well. The shock is seen to 
be spread over several cells (FLASH) or smoothing lengths (GADGET). 
Lower panel as for top panel but including cooling. The analytical solution 
shows that the gas shocks to a lower temperature (due to the smaller differ- 
ence between the incoming gas velocity and the shock velocity), followed 
by a 'cooling tail' in the post-shock region. When simulated using GAD- 
GET SPH particles shock in several steps before reaching their maximum 
temperature. As they do so, particles cool to some extent in the smoothed 
shock and hence reach a lower maximum temperature than the analytical 
solution (the SPH shock is also offset to the left of the analytic shock). In 
the FLASH run gas gets shocked to higher temperatures, closer to the ana- 
lytical solution. Note that as the gas gets compressed the downstream SPH 
smoothing length is smaller than the FLASH cell size. 



we use the higher than usual value of the artificial viscosity (a = 3) 
in GADGET, both the SPH and AMR codes handle this shock well 
(as expected), with the shock smeared out over a few times the reso- 
lution length h in SPH, and a few cells in FLASH. At GADGET'S 
default value for the artificial viscosity (a — 1) we find that this 
is too high a Mach number shock to be handled (we do however 
return to the original value when we study the lower Mach number 
shocks in section[3]). In Fig. [3] we tested both altering the value of 
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Figure 3. Comparison of the effects of altering the viscosity and timestep 
on the shock from the upper panel of Fig. [2] Red lines show the higher 
viscosity (a = 3) scheme with adaptive time steps, green lines the same 
viscosity but a global minimum timestep (set to the minimum Courant step 
of all particles) and purple lines the global minimum timestep but with the 
default viscosity (a = 1), Mack line the analytic solution. Dashed lines 
show the 10th and 90th percentiles. 



the artificial viscosity and adjusting the maximum time step (be- 
tween GADGET'S default adaptive scheme and a global minimum 
Courant step applied to all particles). The higher value of artificial 
viscosity was found to significantly reduce the scatter in the post 
shock thermal energies, superior to a reduction in the global time 
step. We would, however, expect that at very high Mach number 
shocks a more conservative time step would be required. 



25 



20 



15 



10 



_L_ 



1.0 




0.9 




0.8 




0.7 




0.6 
0.5 


< 
< 


0.4 




0.3 




0.2 




0.1 




0.0 





-20 -15 -10 

x/Ax\ 



25 



20 



15 



10 



1 — 

Analytic 

SPH 
AMR 



"T 



QOO OOOOOOOOOOOO^ 1 



1.0 




0.9 




0.8 




0.7 




0.6 
0.5 


< 
< 


0.4 




0.3 




0.2 




0.1 




0.0 





-20 



-15 -10 

x/Ax\ 



2.3.3 Test with cooling 

First let us consider the case of cooling for the M = 4.7 shock. 
This should result in a gas temperature of less than (T± + Tn)/2, 
i.e. we are on the left side of the cooling spike. The initial collision 
of gas can result in higher temperatures and follows an evolution 
for which we have no analytic solution, before settling down to our 
stationary case. 

In Fig[2](lower panel) we see the FLASH and GADGET rep- 
resentations of these shocks. Both codes reach a maximum temper- 
ature which is lower than that of the similarity solution. In SPH we 
attribute this to 'pre-shocking' , i.e particles will shock in several 
stages and cool as they are being shocked. In FLASH we attribute 
this to the cooling operation being applied after the hydrodynamics 
in a time step, such that we do not record the post- shock temper- 
ature. Neither GADGET nor FLASH has the resolution to repro- 
duce the cooling tail particularly well here, although the final cold 
state is achieved in both cases. 

For our second cooling test we look at a more extreme case, 
M = 6.04. This shocks the gas up to a temperature T > T± from 
which it cannot cool, hence the analytic solution is the same as for 
a shock without cooling. In Fig. [4] we show the left hand side of the 
shocked regions. Here the FLASH simulation reproduces the ana- 
lytical result very well, but the GADGET simulation suffers from 
much more severe numerical overcooling through the pre-shock re- 
gion, which prevents the gas from reaching the temperature from 
which it is unable to cool (due to our choice of A(T)). As a re- 
sult we see pile-up of high density cold gas around x = 0, and 



Figure 4. As for Fig.[2]but for an M = 6.04 shock. Upper panel shows 
the case without cooling, with a higher post-shock temperature than Fig. 
[2] Lower panel, the case with cooling. Here the SPH particles shock over 
several smoothing lengths, allowing them time to cool. Unfortunately, this 
means they never reach the higher temperature where cooling vanishes and 
their temperatures decline to the pre-shock value, forming a cold dense re- 
gion similar that in Fig. [2] The shift of the position of the shock front is due 
to the conservation of mass; cooling allows the post-shock gas to be com- 
pressed down to a small region around x = 0. We note here that FLASH 
adequately captures the post-shock temperature even when cooling is in- 
cluded. 

the shocked region is left far behind that of the FLASH rurj^j As a 
result of this overcooling the SPH simulation fails to form any hot 
gas at all. We note, however, that this is a general problem and not 
specific to either GADGET or SPH. 

2.4 Convergence study for GADGET results 

As it stands, we can be confident that the results we have just given 
for SPH are not converged, as they have failed to reach our stable 
analytic solution. The problem we have attempted to solve involves 
no elements that an SPH code would not be expected to handle in 

2 Note that if cooling is prevented, the shock speed will be much higher 
relative to the rest frame. This is easily understood in terms of conservation 
of mass, the gas is shocked to a lower density and a much larger region is 
required to contain it. 
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Figure 5. As for the lower panel of Fig. [2] but for two different SPH particle 
resolutions :red points are the SPH particles as for the lower panel of Fig. 
[2] whereas green points are for a 2x better resolution run (i.e. a factor of 8 
times more particles). The lower resolution run reproduces the temperature 
peak to within 25%, for the higher it is around 15%. When the simulation 
is close to convergence, we would expect AT oc Au oc h, the smoothing 
length, i.e to get within 1% of the temperature we would need a factor of 
~ 10 4 more particles (compared to the lower resolution run). At higher 
resolutions the offset between the exact and simulated shock fronts is also 
reduced. 



the limit where the SPH resolution length h — > 0, and using a good 
prescription for artificial viscosity. The outstanding question here 
is thus only one of how much resolution is required; to this end 
we re-ran the A4 = 4.7 simulation with a factor of 8 increase in 
the particle counQ(from « 80, 000 to « 660, 000). Let us first, 
however, make some general remarks about the problem. 

Given the maximum temperature, T s , of the radiating shock 
we can estimate the error AT of the SPH maximum temperature 
by estimating the radiative cooling over the shock (the physical 
shock is non-adiabatic and so occurs on timescales many orders 
of magnitude shorter than that of the cooling, as discussed in the 
Introduction). Assuming the width of the SPH shock is ~ h, the 
temperature difference will be given, to first order, by 



AT 



voksTo 



■ m p \u A (T s )\ 



h 



Ax A 



(19) 



where we have assumed that all the mechanical energy has been 
converted into thermal energy, and that T s ^> To. For larger 
smoothing lengths we expect AT to become sub-linear in the 
smoothing length, since the cooling is weaker at lower tempera- 
tures (assuming we are on the left hand side of the 'cooling spike'). 

If we apply this argument to Figure [5] we see that increasing 
the number of particles by a factor of ~ 8 (i.e. reducing h by a 
factor of 2) reduces the temperature error by a factor of ~ 1.5, sug- 
gesting that we are not quite in the linear regime. To reach a temper- 
ature within 1% of the analytic temperature would seem to require 
increasing the particle count by a factor of ~ 10 4 . Although this 
is (barely) possible for this particular case, such resolution is not 



3 One can, of course, successively reduce the width of the shock tube in a 
3d simulation to achieve the same scaling as a Id one, i.e. h oc N~Tp U . In 
an astrophysical simulation, however, this option is usually unavailable as 
the shock will be embedded in an environment which needs to be simulated 
in full 3d. 



feasible in cosmological calculations. Most shocks in cosmologi- 
cal simulations will occur at lower resolution than we have used 
in this test case. Therefore we should seek an alternative solution 
involving a switch to prevent cooling during the shock process. We 
intend to explore such a switch in a further paper. 



3 A RESOLUTION CRITERION FOR RADIATIVE 
SHOCKS 

Having established the difficulty of modelling some shock prob- 
lems with radiative cooling we now wish to obtain a criterion 
against which we can judge simulations. Such a tool will allow us 
to identify those simulations where resolution is not a problem and 
those where more care is required. In the following section we will 
discuss the effects of resolution in quite a general way before de- 
riving a metric from a simple model problem. We will frame our 
discussion in terms of SPH, however there will be analogous argu- 
ments for mesh codes. 

Let us take a general case of a numerical simulation of a ra- 
diative shock. We assume we have pre- shock gas with velocity, 
specific internal energy and density v^u^p which passes through 
a shock and comes to rest (v is the velocity of the incoming gas 
with respect to that of the post- shock gas). First we note that the 
SPH shock has a width which is some multiple of the smoothing 
length h oc (mspu/ p) 1 ^ 3 (for a mesh this would be the width of a 
cell), where the numerical factor will include some dependence on 
the artificial viscosity prescription. The change in thermal energy 
will be Au oc v 2 by energy conservation, and thus we can define a 
rate of 'shock-heating', 



Au/ At 

, 3/TOSPHx-l/3 



(20) 
(21) 



where k is some constant depending upon the details of the SPH 
scheme used (e.g neighbour counts). This heating rate is entirely 
numerical, as can be seen by the presence of the SPH particle mass 
mspH: in the continuum approximation of the underlying fluid 
equations the shock heats the gas instantaneously, hence the heat- 
ing rate is singular. As we reduce the particle mass, h decreases 
and the numerical rate at which particles are heated over the shock 
front increases. 

By taking the ratio of the physical rate of gas cooling to the 
numerical rate at which the gas is shock heated (which, ideally, 
we would wish to be almost infinite) we can analyse the effects of 
shock resolution. Only if the absolute value of this dimensionless 
ratio is small do we expect the shock heating to overwhelm the 
cooling, i.e. we require 



/ raspH 



c 3 M 3 \ p 



1/3 



< ?7sph : 



(22) 



for the numerical solution to achieve close to the correct post- 
shock temperature, where 77 is a dimensionless parameter. Here we 
have written the velocity of the incoming gas as v — Mc in terms 
of the Mach number and the upstream sound speed, c = cq. The 
equivalent for a mesh code can be written with the side length h of 
a cubic mesh cell written in terms of the mass enclosed, ttiamr = 
ph 3 , to give 



U A 



C 3 M 3 
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1/3 
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(23) 
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Figure 6. Specific thermal energy vs. position for radiative shocks with a 
Heaviside cooling function, Eq. {24} . Black, blue, green, red crosses are for 
SPH simulations with dimensionless cooling rates of A=0.11, 0.32, 0.53, 
0.74, respectively; positions are quoted in units of the initial (pre-shock) 
particle spacing £±xq, thermal energies in units of the initial thermal en- 
ergy uq . Solid line indicates the analytic instantaneous post shock thermal 
energy u s /uq = 2.44. Dashed line indicates mid-point energy between 
the initial and final thermal energy | (u s + uq)/uq. When the cooling rate 
is low (A = 0.11, black crosses), the numerical overcooling is small and 
the simulation gets close to the right post- shock temperature. Increasing the 
cooling rate degrades the accuracy of the numerical result. We use this to 
set a maximum cooling rate that the simulation can tolerate, for example by 
requiring that the simulated post-shock temperature be larger than half the 
correct result (horizontal dashed line). 



In the subsequent section we attempt to determine a reasonable 
value of r\ which we can use to evaluate other simulations. 



3.1 Heaviside cooling function 

To achieve an extremely simple radiative shock we set up a wall 
shock (see section [2| with a low Mach number M — 2 and the 
piecewise cooling function 



3/2 

Ax 



0, u ^ uo 

-A(^), u <u 



(24) 



where Axo = (mspu/ po) 1 ^ 3 is the initial interparticle spacing, 
A > a dimensionless cooling parameter, uq the initial specific in- 
ternal energy and as in Section[2]we use 7 = 5/3. In the SPH simu- 
lation the particles are initially arranged on a cubic lattice of dimen- 
sions 1024x8x8 in units of Axo (1024 referring to the long x direc- 
tion). The simulations were all performed with periodic boundary 
conditions. Usually a cooling function would be independent of the 
interparticle spacing, however we chose to re-use our initial condi- 
tions whilst adjusting the dimensionless constant A, and in this way 
scale the LHS of Eq.|22] This is equivalent to using a fixed cooling 
function but adjusting the interparticle spacing. 

We now make a couple of observations. Firstly, we note that 
we can calculate the instantaneous post- shock state using the equa- 
tions derived in Appendix Al to find p s /po 



2.52, ; Us/uq — 
2.44. We note that this ratio is independent of the cooling parameter 
A. Increasing A in the simulations, however, we expect the maxi- 



mum post-shock temperature to fall as thermal energy is radiated 
away over the numerically broadened shoclQ 

In Fig. [6] we see the results of these simulations plotted at a 

l— ' —1/2 

time of t = 141?x Axo. We note that the particle distribution 
in the pre-shock region has also diverged from a lattice arrange- 
ment (if it were still a lattice the particles would appear at multiples 
of Axo) into something more glass-like. This is to be expected as 
the SPH equations of motion favour a large distance to the nearest 
neighbour for a given density, which can be acheived with a close- 
packed or glass-like arrangement. The position, velocity and Mach 
number of the shock at late times are independent of the cooling 
function, for fixed uq (provided the cooling function restores the 
thermal energy of the gas to uq) as is shown in Appendix [Al] 

With a low cooling parameter (A = 0.11, black crosses) we 
see that the post-shock thermal energy reaches near the theoreti- 
cal value, whilst with a high cooling parameter (A = 0.74, red 
crosses) we see that the simulation produces almost no hot gas. We 
take the mid-point of the thermal energies as a minimum value the 
code should reach to give at least approximately the correct answer. 
From Fis. [i]this corresponds to a cooling parameter of A « 0.4. 
Substituting this maximum value into Eq. ( |22| allows us to evaluate 
the parameter 77 as 



3/2 



1 



c 3 M 3 



2/3 



A(7(7-1))" 3/2 A^" 
0.08, 



2/3 



(25) 

(26) 
(27) 



(where we have used the post- shock density p — p s — 2.52pn), or, 
using Eq. ( [22] ), 

1/3 



\U A \ 



( msPH 



< 0.08. 



(28) 



c 3 M 3 V p 

This is the value of 77 we will use throughout the remainder of this 
paper. 

A similar analysis with FLASH yields an only slightly weaker 
criterion, 

1/3 



1 ( m A MR A 

p ) 



c 3 M 



< 0.09, 



(29) 



where ttiamr refers to the mass contained in a mesh cell (since the 
mass in cells varies we have taken ttiamr to be the value in the 
cell immediately to the right of the shock, for consistency with the 
evaluation of p). 

It is worth discussing the differences between a grid and an 
SPH scheme when the adaptive capabilities are utilised. SPH has 
a resolution (smoothing) length which refines in areas of high den- 
sity as p -1 ^ 3 . AMR on the other hand can have much more general 
refinement criteria, for example allowing higher resolution to be 
applied on features which need not be dense (e.g. shocks). As such 
AMR has something of an advantage when it comes to shocks, as 
almost all refinement schemes will refine over discontinuous vari- 
able to the maximum level, and hence there is no need to impose 
the refinement criterion Eq. {29). Of course it is possible that a re- 
gion of space in the simulation is already maximally refined, yet 

4 One might have expected the post shock thermal energy ratio for a mach 
2 shock to be precisely 27(7 — 1) + 1 = 3.22 as in the case without 
cooling, however the immediate post-shock region is still in motion w.r.t 
the final cold post-shock gas (hence A = is a special case). 
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Figure 7. Cooling functions used in the GIMIC simulations at redshift 0. 
Upper solid and dot-dashed lines represent astrophysical cooling functions 
for a plasma with metallicity [Z/Zq] = 0,-3 (where square brackets 
denote the base- 10 logarithm) respectively in the presence of an ionizing 
background (Wiersma et al. 2009a I. Lower solid line shows a cooling spike 
such as in section[2] for comparison, on the same logarithmic scale. Dotted 
line shows cooling due to oxygen only, again assuming a solar abundance. 



even so fails to satisfy the criterion Eq. ( |29| >. We can then interpret 
this as a test of how well the finite resolution of an AMR simulation 
represents the physics in the problem. 

To refine a simulation in a given volume V of a 2 dimensional 
structure (such as a shock) to scale h in SPH requires iVspH oc h~ 3 
particles, whilst in AMR one would only need N ce \\ oc h~ 2 cells 
(we note that limitations on the refinement level between cells does 
not in general alter this scaling relation). 

This can be contrasted with a sheet like structure in a vacuum 
(e.g. a gravitationally collapsed disk of thickness <C h), which will 
only require A^sph oc h~ 2 particles, the same relation as AMB0 

Note that we have been concentrating on how well the simu- 
lations reproduce shocks in the presence of cooling. In practise we 
would also like the code to correctly reproduce the cooling tail, i.e. 
the cooling of the gas once it has passed through the shock (the 
right-hand side of Fig. [6). Clearly here we would like to resolve the 
cooling length from Eq. |l7| , by requiring that Ax a ~ h in SPH 
(or the cell size in mesh codes). 

In the subsequent section we will apply the resolution criterion 
we derived to estimate in which areas of cosmological simulations 
numerical overcooling could be problematic. 



4 EFFECTS OF RESOLUTION ON GALAXY 
FORMATION 

4.1 Galaxy formation simulations 

In this section we apply the resolution criterion of Eq. (28} to differ- 
ent regions of temperature and density in the the GIMIC SPH sim- 
ulation of Crain et aL] (2009) . First we will plot the distribution of 
gas in temperature-density space and identify some environments 
of interest. We will then discuss the radiative shock resolution in 
this parameter space, but also explore mitigating factors which may 

5 As such SPH could be viewed as a refinement scheme optimised for grav- 
itationally collapsed structures. 
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Figure 8. Resolution requirements for correctly representing shocks in dif- 
ferent regions of a temperature-density diagram. Solid contours are labelled 
with the minimum values of the shock speed, v = cA4, obtained from 
Eq. {28) , that avoids excessive overcooling in the shock precursor, for sim- 
ulations using an SPH particle mass of raspH = 1O 6 M0. The radiative 
cooling rate adopted is that of a plasma with solar abundances of elements, 
[Z/Zq] = 0, as shown in Fig. [7] The overlaid red shaded region is the 
phase density in (T, tih) space of SPH particles in a cosmological feedback 
simulation (see text). Bold letters refer to example environments described 



in Table ^ Heavy black and grey dashed lines refer to t coo i 
t C ooi — respectively. See text for discussion. 
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allow us to have confidence in simulations even when they fail to 
accurately resolve the shocks. 

The GIMIC simulations are zoomed re- simulations of nearly 
spherical regions picked from the Millennium simulation ( Springel 
|et al.|2003] >, including gas dynamics. Each sphere has a radius of 
ISh 1 Mpc, and the SPH particle mass is m S p H « 1O 6 /i _1 M . 
Radiative cooling in the simulation includes line cooling of eleven 
elements, Compton cooling with the CMB and thermal brem- 
strahlung in the presence of a uniform but evolving ionising back- 
ground, as described in Wiersma et al. 2009a (see Fig. [7]). The back- 
ground cosmology, as for the Millennium simulation, is Q m — 
0.25, Q A = 0.75, Q h = 0.045, n s = 1, a 8 = 0.9, H Q = 
100 h km s _1 Mpc -1 , h = 0.73. The enrichment of gas by nu- 
cleosynthesis in stars is described in Wiersma et al. (2009b). Photo- 
heating, radiative and adiabatic cooling, shocks induced by galac- 
tic winds and due to accretion, result in gas occurring over a wide 
range of densities and temperatures, illustrated in Fig. [8] Five points 
A-E in this T—p space correspond to physical states where we want 
to investigate to what extent the simulation properly resolves radia- 
tive shocks (see also Table [TJ. For a general discussion of these 
diagrams see Wiersma et al. ( 2010). The simulation code described 
here was also used in the OWLS project ( Sch aye et al.|2010| ). 

Point A is a typical IGM point outside virialised halos, at low 
density and temperature. Here we see a very well defined mild up- 
ward slope of temperature with p of the post reionization gas. This 
gas is cooling due to adiabatic expansion of the universe and is 
being photo-heated by the UV-background . For a recombination 
coefficient oc T -0 ' 7 this will at late times result in a temperature 



1/17 |Hui&Gnedin|l997 



Theuns et al. 



density relation of T 
[T998] ). 

Point B corresponds to gas heated in an accretion shock, 
falling into a galactic halo, or shocked by a galactic wind. Mechan- 
ical energy has been converted into thermal energy, and the density 
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Table 1. Astrophysical shock environments identified in Fig. [8] 



will jump by a factor up to ~ 4. When this gas cools, it will form 
the warm gas of point C which may condense to fuel star formation 
in a central galaxy ( White & Rees 1978 ). 

On the far right, nn > 10 _i cm -3 , is a sharp vertical feature 
in the distribution of SPH particles. This boundary delineates cold 
halo gas from gas which undergoes star formation in the model 
used in GIMIC. The denser gas represents a multi-phase inter- 
stellar medium, for which the imposed pressure-density relation in 
GIMIC is p oc rij 1 / 3 , known as an effective equation of state for 
the ISM. Such a state is intended to mimic the physical pressure 
response in dense gas undergoing star formation (point D), com- 
pressing this gas results in significant star formation with associ- 
ated feedback (see Schaye & Dalla Vecchia (2008 ) for motivation 
and details). The SPH density then represents a volume average 
density of star-forming gas, whereas the physical ISM lies in ap- 
proximate pressure equilibrium, but with a hot and cold phase and 
corresponding variation in densities. In particular the simulation 
does not allow this gas to cool radiatively. Finally, point E repre- 
sents the domain of type II SNe that ignite in the hot (10 6 K) sparse 
phase of the ISM, generated by the activity of previous generations 
of SNe. We note that there is little gas marked in this phase as the 
cooling time is short. 

Now let us consider this simulation in terms of its ability to 
resolve the radiative shocks that occur in these 5 regions. Equation 
(28} suggests that a radiative shock of velocity v will be resolved if 
it satisfies 



v > 



\UA[ 

0.08 



1/3 / \ 1/9 
/ msF h x 

V P 



(30) 



i.e. there is a minimum shock velocity which can be resolved. 
Shocks below this velocity will tend to artificially radiate away 
their energy because there will be cooling through the (artificially 
extended) shock region. Shocks above this velocity will heat up the 
gas on such a short timescale in the simulation that the cooling in 
the shock region will make little difference to the final result. 

In Fig. [8] we plot contours of given v, the minimum shock ve- 
locity for which there is no significant overcooling in shocks. At 
each temperature and density a cooling rate is evaluated (using the 
cooling rate from Wiersma et al. 2009a shown in Fig. [7] evaluated 
at solar metallicity, and in the absence of an ionising background), 
which is combined with a particle mass of ttisph = 10 6 M©, to 
derive a minimum shock velocity which can be resolved. Note that 
Eq. (30} is very weakly dependent on particle mass, and thus chang- 
ing mass resolution is a very ineffective way of shifting the con- 
tours. These contours represent the minimum velocity shock which 
can be resolved at each T, p. Any shocks at lower velocities will 
appear artificially colder due to resolution effects. 

A key point, however, is that even if we fail to resolve the 
radiative shock, the cooling of the gas in many cases may be in- 
evitable anyway. Indeed, there can be situations where other pro- 



cesses are occurring on much longer timescales than the cool- 
ing, and for which having an incorrect thermal history of the gas 
is not a problem as far as the dynamics of the system is con- 
cerned Establishing a general criterion for these is not trivial, 
here we will simply compare to the locally estimated dynamical 
time tdyn = (Gp)~ x ^ 2 as indicative of the timescales for other 
processes. We assume that the simulation will cool adequately if 
tdyn ^> tcooi even in the case that radiative shocks are resolved 
poorly (we define t coo i = \ua \ /u). The heavy dashed black con- 
tour in Fig. |8]represents the line where tdyn — tcooi • All points to 
the left of this represent tdyn < t coo i so certainly we would wish 
to completely resolve any shocks here. We have also included in 
dashed grey a line where 0.1 t^n = t coo i to demonstrate a some- 
what stronger limit. The necessity of resolving the thermal history 
to the right of this line, is questionable, because the gas cooling 
time is so small in any case. Of course these simulations assume 
ionisation equilibrium and optically thin gas, so the cooling rates 
may have been overestimated. In addition if one were to attempt 
to track the chemistry of the shocked gas, for example the forma- 
tion and destruction of molecular hydrogen, then having a correct 
thermal history could still be important ( |Abel et aLfT 997 ). 

Now let us evaluate the resolution criterion of Eq. |30|) for the 
five diverse environments of Table [T] First we take point A, for 
radiative shocks in the IGM; here we can see that the low density 
and cooling rate combine to allow us to resolve all shocks above 
5 km s _1 , almost certainly much exceeding our requirements. 

For point B we have taken a value for gas heated by a virial 
shock to 2 x 10 6 K. The higher density and cooling rate here 
push up the minimum shock velocity we can resolve to around 
100 km s _1 , comparable to the virial shock velocities ^200 them- 



^200 = (10GH(z)M 200 ) 
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( |Mo e t al. 1998), where G is the gravitational constant, H(z) the 
Hubble parameter and M 2 oo is the virial mass of the halo. For 
z > 1 we can approximate the Hubble parameter as H(z) « 
#o^m 2 (l + z) 3/2 to see 
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(34) 



where Q m is the matter density in units of the critical density and 
h = i/o/100 km s _1 Mpc -1 . For lower mass haloes the gas actu- 
ally cools even faster and the shocks are more difficult to resolve, 
however the cooling may be so short as to save the situation. We 
explore this situation further in Section [43] 

For point C we consider the warm gas within galactic disks. 
The minimum shock velocity which can be resolved close to the 



6 Note that even if dynamics is not affected, there may be other conse- 
quences of numerically underestimating the amount of hot gas, for example 
when calculating the spectrum of cooling radiation. 
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Figure 9. Contours of maximum SPH particle mass raspH required to pre- 
vent numerical overcooling at a virial shock, for gas ([Z/Z@] = — 3) ac- 
creting onto haloes of different virial masses M200 at a given redshift z. 
Coloured lines corresponding to raspH = 10 s , 10 6 , 10 4 , 1O 2 M lim- 
its are represented by the thin maroon, yellow, cyan and blue lines respec- 
tively. Black lines compare cooling time to the dynamical time of the halo: 
tcool = tdyn (heavy solid line), t coo i = 2t dyn , t cooi = \t dyn (heavy 
dashed lines). The shaded grey region denotes t coo i > tdyn- Numeri- 
cal overcooling due to lack of resolution in regions where t CO oi ~ tdyn 
will likely affect the dynamics of the accreting gas, hence SPH simulations 
would appear to need resolutions ~ 10 6 M@. 



star formation threshold is higher than for point B, because the 
cooling rate is higher. However the cooling time is so much smaller 
than the dynamical time in the disk, so that we suspect that gas 
cooling will be inevitable in any case. This suggests that, although 
numerical overcooling is potentially severe here, it is unlikely to 
have any important effects on the evolution of the disc. 

At a higher temperature than C we have point E, a fidu- 
cial point for a (type II) supernova blast wave shocking to hot 
(~ 10 6 K), rarefied (n H ~ 10" 2 cm" 3 ) ISM (irradiated and in- 
flated by the massive progenitor star). We have a high minimum 
resolved shock velocity due to the fast cooling of this gas, making 
its simulation problematic. We expect the gas to form a cold, dense 
shell (Co x|1972) , and the lack of resolution to manifest itself pri- 
marily in an alteration of the onset of this phase. We discuss the 
implications for SN feedback further in Section [43] 

Finally, for point D we have considered the case of an AGN 
outflow shocking into a dense ISM of uh ~ 1 cm -3 . The mini- 
mum resolved velocity is so high here that we can have little con- 
fidence in the simulated shock dynamics (excluding the most basic 
properties such as conservation of momentum). The gas is cooling 
fast compared to dynamical time scales, yet we have similar con- 
cerns to point E about the artificial suppression of feedback. 

4.2 Virial Shocks 

We now consider the effects of the resolution requirement Eq. (28} 
on the discussion of cold accretion and virial shocks around haloes. 
Here we are following the ideas of spherical collapse set out in 
White & Rees (1978 ). The basic question here is what is the fate 
of gas as it accretes on to a halo, and gets shocked as it converts its 
mechanical energy into thermal energy. If the cooling time is short, 
then this hot phase will be a temporary one, however if the cooling 
time is long, then a hydrostatic hot halo of gas will form within the 
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Figure 10. As for Figure [9] this time including a uniform ionising back- 
ground, see text for details. In the gold region (lower left of the figure), the 
gas is being heated rather than cooled, so resolution of the shock is of lesser 
importance. 



halo, the scaling relations for which can be found in e.g. |Mo et al.| 
(1995) . 

The properties and stability of such spherical shocks have 
been studied by e.g. Birnboim & Dekel (2003 ), depend on mass 
and redshift. More massive haloes have hotter shocks with longer 
cooling times. At a given mass lower redshifts imply lower densi- 
ties and hence slower cooling, and hence easier build-up of a hot 
halo. It must be recalled, however, that in this situation geometry 
will also play a role. If the gas accretion can achieve a configura- 
tion where it will penetrate farther into the halo (e.g. in filaments), 
it will shock at higher densities and generally have a shorter cooling 
time (the situation is complicated by the fact that the gas continues 
to accelerate, and so the shock will generally be hotter). 

Here we first consider applying our resolution criteria to the 
spherical case. Assuming a spherical halo of mean density 



P200 = 200 



(3H 2 \ 
\8tvGJ 



(35) 



and virial mass M200, we take the accreting gas to shock to the 
virial temperatur^T2oo and virial velocity V200 given in Eqs. 
and (33}, respectively. For the baryon density at the edge of the halo 
we use 

1 Q_ b 

3 n 



Pb = ~- 



"P200, 



(36) 



where the factor 1/3 is the ratio of edge to mean densities for an 
isothermal sphere of profile p — po(r /ro)~ 2 . We can then apply 
our shock resolution criteria in terms of the maximum mass of SPH 
particles that do not suffer from numerical overcooling in the shock, 



m S PH,max = ^ 3 |^A(T 20 o)l ^200 Pb , 



(37) 



where our convergence tests suggest that rj « 0.08. 

Equation ( [37) defines curves in a plot of virial mass M200 
versus redshift z, shown in Fig. [9] for a cooling rate appropriate 
for a plasma with solar abundance ratios but mean metallicity of 



7 The infalling gas has actually twice this energy so if it shocks into the 
rest frame of the halo the temperature will be increased by a factor of 2, 



however we will ignore this for now. 
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[Z/Zq] — —3 (we have chosen the lower metallicity as more 
representative of accreting gas that has yet to be enriched by sev- 
eral generations of star formation). In Fig. [lO]we show the case 
where cooling is partly suppressed by the presence of a uniform 
ionising background (see Wier sma et al.|p00"9a| ) for details). Each 
thin coloured line represents the Umiting particle mass required to 
prevent numerical overcooling in the corresponding spherical ac- 
cretion shock. Clearly the resolution requirement is punitatively 
strict (smallest rasPH) for small haloes (M 20 o ~ 1O 8_1O M ), 
especially at high redshifts (z ~ 9). Intuitively this can be un- 
derstood because these masses correspond to virial temperatures 
near the peak of the cooling function, and at higher redshifts the 
mean baryon density (and thus collisional cooling rate) grows. In 
the presence of an ionising background cooling is suppressed in 
lower-mass haloes that cool mostly through hydrogen lines. At 
lower masses the ionising radiation has a very large effect, and gas 
will be photo-heated instead of cooling (Oka moto et al.|2 008). 

However, even though lack of resolution will lead to overcool- 
ing in some haloes, the cooling time in these haloes may be so short 
that the gas would cool quickly anyway. The grey area in the figure 
indicates where the dynamical time in the halo is shorter than the 
cooling time: in this region we expect that numerical overcooling 
may prevent the formation of a hot halo. Conversely, in the white 
region, cooling is so fast that even if a hot halo were to form, it 
would quickly cool. The demarcation line between these scenarios 
follows closely the ~ 1O 6 M0 limiting SPH mass (yellow line). 
Simulations run with that resolution or better will be able to form 
hot spherical haloes in situations where we would expect such a 
hot halo to form. At lower resolution, simulations may artificially 
suppress the formation of a hot halo due to numerical overcooling 
in the accretion shock. 

Our considerations apply only at the virial radius. However, 
nearer the centre of halos we expect this conclusion to remain valid, 
as the cooling time diminishes faster than the dynamical time. The 
very high mass halos have 2td yn < t coo i (heavy dashed black line) 
and we expect these to be in near hydrostatic equilibrium. As a re- 
sult of these analyses we conclude that 1O 6 M is a sensible upper 
limit for the gas particle mass in cosmological simulations intend- 
ing to capture the evolution of proto-galactic haloes, although lower 
masses enable more accurate resolution of the thermal history of 
gas in lower mass haloes. 

4.3 Thermal feedback 

Thermal feedback refers to the mechanism whereby injection of 
thermal energy into the ISM causes adiabatic expansion of the gas 
and subsequent suppression of star formation due to the diminished 
density. The simplest model to envisage is perhaps that of a single 
supernova creating a hot, spherical, rarefied, blast wave. On larger 
scales, however, we expect to see analogous effects from star form- 
ing regions and AGN. In this section we intend to consider our 
results in terms of thermal feedback in SPH. We will review the ba- 
sic physics of thermal feedback and its role in galaxy evolution. We 
will then discuss its implementation in SPH and derive some quan- 
titative criteria for accurately resolving it. Finally, we will relate 
our observations to the feedback experiments in other work. 

We begin by considering the problem of simulating a super- 
nova blast wave. Here we are primarily concerned with the situa- 
tion where we may artificially radiate away the thermal energy of 
the blast wave due to a lack of resolution. This would result in the 
premature transition from a thermally driven to a momentum driven 
phase. 



A concise overview of the evolution of a supernova remnants 
can be found in Cox (1972 ). Essentially the blast wave will follow 
a Sedov-Taylor self-similar solution ( |Sedov|1 959 ) until the shock 
temperature T s falls to a value where the radiative cooling exceeds 
the cooling via adiabatic expansion. A full calculation is beyond 
the scope of the present paper, however, we can get close just by 
dimensional considerations 



k B T s = (Eimln A n k Q ) ini 
« fc B -4xl0 6 K, 



(38) 
(39) 



where E — 10 erg is the S Ne energy, nn = 1.0 cm , A = 
10 -22 erg cm 3 /s. The value of |Cox|l972| is a factor 2 smaller, at 
2.0 x 10 6 K. 

In a simulation we would like to resolve the transition in the 
supernova shock from being pressure driven to being momentum 
driven, which typically occurs for shock temperatures T s ~ 2.0 x 
10 6 K; numerical overcooling may cause the shock to transition too 
early, hence underestimating the feedback effect of the explosion 
on star formation in the surroundings. Using the Sedov similarity 
solution for a 3 dimensional blast wave in a uniform cold medium 
of adiabatic index 7 = 5/3 and density po (which we will shortly 
take to be the density of the ISM), we can then write the pressure 
and temperature just inside the shock wave in terms of the shock 
velocity v s , as 



Ps 

k B T s 



7 + 1 
7 - 1 

2 2 
7 + 1 

-Pa 



o 7-1 - 2 
2 - — — —fJ>v 8 , 



(40) 
(41) 
(42) 



Ps (7 + I) 2 ' 

where ft is the mean particle mass. Combining this with the resolu- 
tion criterion of Eq. ( [28] ) we find that (excluding fairly pathological 
cooling functions) the cooling will be hardest to resolve at the lower 
temperatures, i.e. atT s = 2.0 x 10 6 K. 

Applying fiducial values for the ISM of fx « 0.6 m p at T s = 
2.0 x 10 6 K yields a shock velocity, density and cooling rate of 

Vs = 380 km s" 1 (43) 
p s « 9 x 10" 24 gem" 3 (44) 
u\ A « -180 cm 2 s" 3 , (45) 



and the corresponding radius of the blast wave 



1.15 



Po 



1/5 



2/5 



1.15 5/3 
15 pc , 



Po 



2/3 



-2/3 



(46) 

(47) 
(48) 



where E ~ 10 51 erg is the thermal energy injected by the explo- 
sion. The corresponding limiting SPH particle mass that avoids nu- 
merical overcooling, evaluated from Eq. ([28]), is 



mspH = 70M G 



9 x 10- 24 gcm- 



( v -± y 

V380 kms" 1 / 



U A 



180 cm 2 s- 



(49) 



The small values of both the radius and the required mini- 
mal mass resolution imply that most cosmological simulations of 
galaxy formation are far from resolving individual SN explosions, 
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however detailed simulations of high-z dwarf galaxies do indeed 
already reach such extreme resolutions {e.g. Wise & Abel (2008)). 
For a state of the art mass resolution for a cosmological simulation 
of say 10 5 M©, a star particle really represent very many stars and 
hence also many SNe. Simply scaling-up Eq to represent the many 
SNe that go off does not really help much, as for example the blast 
radius only scales oc E 1 ^ 3 . In reality different SNe will go off in 
different places, and once the density of the ISM is decreased due 
to one explosion, another explosion in the lower density gas will 
have a much larger effect, eventually resulting in a percolating hot 
phase. 

However these small scales cannot yet be resolved in current 
cosmological simulations, hence they fail to follow the transition 
from pressure to momentum-driven SN shells. One can try to model 
the expected effects by simply heating a small number of neigh- 
bouring SPH particles. In this case our resolution study indicates 
that the reheating temperature must be sufficiently high for numer- 
ical overcooling not to affect the dynamics. We can thus generalise 
the above calculation to find the minimum temperature for resolved 
thermal feedback for a given SPH particle mass. To perform this 
calculation we will need to associate a shock velocity with a single 
particle, which we will take to be the Sedov shock velocity for a 
blast wave of the same mass as the SPH particle 



1.15 



5/2 



37(7 ■ 



l.lc s . 



Combined with Eq. \30\ we find 



-2/3 



l.l 2 - 7(7-1) V P 



/ msPH \* ,..2/3 
«A 



(50) 
(51) 

(52) 
(53) 



the minimum thermal energy required to avoid numerical overcool- 
ing. 

Dalla Vecchia & Schaye (2008) argue that for thermal feed- 
back to be effective requires that the sound crossing time across 
an SPH particle, t s = h/c s , be smaller than the cooling time, 
t c = u/\ua\- Using ph 3 ~ ttisph, t s < r c requires that 



^fb • 



raspH 



2/9 



2/3 



(54) 



which is identical apart from a numerical factor to Eq. ( [53] ). Our 
criterion is stronger as it takes into account that the code will in 
practise overestimate the cooling of the gas in shocks; the lower 
value simply requires there to be a shock. 

In Fig.[TT]we explore the parameter space for modelling ther- 
mal feedback in an SPH simulation with 1O 6 M particles. At each 
density, there is a minimum temperature required to drive an adi- 
abatic blast wave phase. The light grey region is defined by the 
sound crossing time argument of Eq. (54) and the dark grey is 
from Eq.([53]>. In the white region we expect effective thermal feed- 
back; in the dark grey region it will be suppressed by overcooling 
in shocks, and finally in the light grey the code will be unable to 
produce a shock at all. 

It is helpful to introduce some numbers. For gas with hydro- 
gen number density nn = 1 cm -3 , a raspH = 1O 6 M particle 
would need to be heated to a temperature of T s « 5 x 10 6 K to 
be in the pressure driven phase, according to Eq. (54} . However 
our resolution study suggests that we need a higher temperature of 
« 10 7 K to prevent excessive overcooling through the shock, imply - 




-4 -3 -2 
log 10 n H (cm- 3 ) 

Figure 11. Minimum re-heating temperature T required to avoid numerical 
overcooling as a function of hydrogen number density tih, assuming an 
SPH resolution of raspH — 10 6 Mq and solar metallicities. Cooling is 
so rapid in the shaded regions that the transition from thermally driven to 
momentum driven expansion phases of supernova bubbles is so fast that 
much of the injected energy will be lost to radiation. The light grey shaded 
region corresponds to Eq. {54) , the dark grey shaded is the more demanding 
Eq. ( |53) . The white region is where the reheating temperature is sufficient 
to force thermal feedback despite resolution concerns. Dashed line is an 
estimate of the specific energy of SNe from Ka y et al.|(2002) . 



ing that at this resolution the simulation cannot properly represent 
the effects of thermal feedback. Note also from Eq. ( [53] ), that this 

2 /9 

improves only very slowly with improved resolution, oc m S p H . Re- 
lating back to models of feedback, this is somewhat problematic as 
the specific energy of supernovaejjis estimated to only be around 



2 x 10 K (dashed line in Fig. [ill which can still be too low for 
the simulation code to properly follow the thermal evolution of the 
explosion. 

Clearly this has important consequences for prescriptions 
for supernova-driven thermal feedback. In densities above nn = 

1 cm -3 the thermal feedback starts to reduce its effectiveness, even 
if we ignite all our SNe in a single timestep, yet this is one of the 
key environments where feedback is required. 

One way to drive feedback at this resolution, whilst still main- 
taining a globally consistent initial mass function, is to stochas- 
tically inject the energy of SNe due to star formation. Since the 
mean thermal energy of a single particle can then be greater than 

2 x 10 7 K we can remain in the effective region of Fig. 11 To bal- 
ance the IMF other star forming particles will need to receive less 
or no SNe energy. The alternative is to increase the resolution, but 
again the low exponent of raspH in Eq. ( [53] ) makes this quite pro- 
hibitive. 

The issue is further complicated by the existence of a multi- 
phase ISM not resolved by the simulation. This is the motivation 
for many of the prescriptions for feedback, such as applying a frac- 
tion of the supernovae energy as kinetic energy ( Navarro & White 
1993), disabling the cooling of thermal bubbles ( Gerritse n|1997| > or 
releasing the energy of many accumulated supernovae in one step. 
A more thorough discussion of all these methods can be found in 



8 The specific energy of supernovae is the energy released by supernovae 
per unit mass of star-forming gas ( Kay et al. 2002), i.e. if all supernovae 
were to ignite simultaneously we would reach this mean temperature. 
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Kay et al. (2002). Another approach is to model the net feedback 
effects by a subgrid model (for example an imposed equation of 
state without cooling as in Schaye et al. 2010 ), or to model the hot 
and cold phases separately by representing clouds in the cold phase 
as collisionless particles pooth et al.|2007| ). 

At lower densities it becomes easier to thermally drive a blast 
wave due to the reduced cooling rate, reinforcing the importance 
of simulating a multiphase ISM. For their star formation thresh- 
old of 77-h = 10 _1 cm -3 , the high-resolution OWLS simulations 
of |Schaye et al.| ( [2"010| ) can represent thermally driven SNe at the 
edges of discs, but not in more central regions. Indeed, at higher 
densities the required temperatures rapidly reach extreme values. 
|Booth & Schaye| ( [2009| ) note that in their simulations, AGN feed- 
back requires reheating temperatures T s > 10 8 K, as at lower tem- 
perature the energy is simply radiated away. We believe that this 
problem is not a physical one but one of resolution. At lower tem- 
peratures the density is higher, the cooling faster, and the cooling 
region behind the shock cannot be resolved, as is clear from Fig.[TT| 



4.4 Shocks at the sound speed 

As an interesting aside it is worth considering that there will usually 
be an upper limit to the resolution required. If we assume that the 
weakest shock has a velocity on the order of the sound speed, c s = 
(7(7 — ^) u ) 1 ^ 2 1 men me minimum requirement for the particle 
mass for a given problem will be 



mspH 



(0.08) 3 ( 7 ( 7 -l)) 9/2 min{|^ A |-V /2 p} . 



(55) 



Unfortunately such a limit will usually be very small indeed, at 
least for cosmological simulations, because of the low sound speed 
of cold, dense gas present in galactic discs. However if one chooses 
to go down this path, then one can examine the following criteria. 
If we have a conventional collisional cooling function then ua/p 
is independent of density, and we can make the additional assump- 
tions that ua —> as u — )> 0, and for large u 



\UA[ 
P 



1/2 



(thermal bremsstrahlung), giving 

{p- 2 } , 



mspH oc mm 



(56) 



(57) 



i.e. the smallest particle mass is determined by the highest density 
in the problem. This analysis is of course not valid with Compton 
cooling via the CMB, or the presence of a UV background, as nei- 
ther process is collisional. 



5 CONCLUSIONS 

In this paper we have examined the role of radiative cooling in 
shocks. We have found a general analytical solution for Id piece- 
wise linear collisional cooling functions and compared it to numer- 
ical simulations of the same shock, performed with an SPH code 
(Gadget) and an AMR code (FLASH). These codes smear out 
the shock over several particles or cells, and such an artificial 'pre- 
shock' results in numerical overcooling which may prevent the for- 
mation of a hot post-shock region. We have estimated a general 
resolution criterion to avoid such overcooling, and applied it to the 
problems of virial shocks and the production of hot gaseous haloes. 
We have found that to avoid numerical overcooling of accretion 



shocks onto haloes that should develop a hot corona requires a par- 
ticle or cell mass resolution of 1O 6 M0 (Fig.|lO|), which is within 
reach of current state-of-the-art simulations. 

Similarly, we have applied our estimates to thermal feedback 
from AGN or supernovae blast waves, in the presence of radiative 
cooling. We have seen that the energy required to drive thermal 
feedback at a given mass scale, for current numerical results, is an 
order of magnitude higher than one would expect just from physi- 
cal considerations. For cosmological simulations (1O 6 M gas par- 
ticles) of an 71h = 1.0 cm -3 interstellar medium we see (Fig. flip 
that temperatures in excess of 10 7 K are required to effectively 
drive thermal feedback by avoiding spurious suppression of the 
feedback by numerical overcooling. 

Although all of these issues can be rectified by increasing the 
resolution, the minimum thermal energy of injected feedback re- 
quired to avoid artificial coolin g sca les weakly with decreasing par- 

Consequentially, a potentially 



tide mass, oc mg(, 9 H , see Eq. (53 



fertile region of study may be that of cooling switches, i.e. a cri- 
terion for disabling cooling through a shock. Such a switch would 
allow a simulation to resolve temperatures much closer to the phys- 
ical temperatures of radiative shocks without requiring extreme res- 
olutions. Unfortunately it is not a straightforward problem to have 
a criterion that will consistently suppress cooling in the presence 
of shocks yet does not affect cooling in regions where there are no 
shocks. Since we can never hope to completely remove resolution 
effects it seems sensible to have a more limited aim, perhaps to cap- 
ture the temperatures of shocks up to some maximum cooling rate. 
As such one might wish to suppress cooling, when the cooling time 
is greater than some fraction of the shock heating time. We intend 
to explore this avenue in a further paper. 

Further work could include the effects of shock-induced non- 
collisional ionizational equilibrium (CIE) or non-thermalised gas. 
Since the resolution can make such a significant modification to the 
thermal history of a gas, we expect a criterion due to non-CIE may 
be quite strict. 
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APPENDIX A: RADIATIVE SHOCKS WITH PIECEWISE 
LINEAR COOLING FUNCTIONS 

Al Similarity solution for a Id radiatively cooling shock 

We start with an ideal gas with adiabatic index 7 

P = (7 - l)f>u , (Al) 
and a collisional radiative cooling function 

du\ A = -pf(u)dt. (A2) 
These combine to give an evolution of 

du = (7 - l)-dp - pf(u)dt . (A3) 

Stationary solutions of a post shock cooling region satisfy in- 
tegrals of the mass and momentum equations, i.e. 



p(v-u s ) = po(v -u s ) 

P + p(v-U s ) 2 = po + po(v - U s ) 2 , 



(A4) 
(A5) 



where u s is the shock velocity and p , p , v denote the density, 
pressure and velocity at some arbitrary downstream point. Thus the 
density, velocity and thermal energy can be written in terms of a 
similarity variable A 



p/po 

V - U s 
Vo — U s 

u/uo 
with 

_ Po(v 



(a+l)A _1 -a\- 



Po 



(A6) 
(A7) 
(A8) 

(A9) 



Now we assume we have a piecewise linear cooling function, 
i.e. we can solve each segment separately with the linear cooling 
function 



f(u) — A(u - u c ) . 



(A10) 



where A is some constant and u c denotes the 'cold' thermal energy 
where cooling vanishes. This gives an o.d.e for x of the form 
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which can be solved generally, however in the case of u c 
have the particularly simple case, 
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the left hand limit for A coming from entropy considerations. An 
example cooling shock of this form can be seen in Fig.[T] 



A2 Colliding gas 

Assume two homogeneous flows collide from the left and right, 
with properties po,Po, ±^o- With no cooling a hot, static region is 
created in the centre, with properties p c , and p c . The mass, momen- 
tum and energy equations are 

(vo - u s )po = -p c u s (A12) 
(vq - u s ) 2 pQ + po = p c u 2 s + p c (A13) 

^ + i( 7 -l), 2 = (A14) 
po 2 p c 

where u s is the velocity of the left moving shock in the rest frame. 
Eliminating p c , p c gives 

u s + \ (7 - S)u s v = ^ + \ (7 - 1)^0 , (A15) 

z po Z 

so 



«. = -|(7 - 3)«o - \^vl ( 7 + l) 2 + 16^ . (A16) 

Assume now that there is cooling and that the gas in the cen- 
tre cools to the temperature of the pre- shock gas (where cooling is 
assumed to vanish). In this case the mass, momentum and energy 
equations are 

(vo - u s )p = -p c u s (A17) 
(vq - u s ) 2 po +po = p c u 2 s + p c (A18) 
Po/po = Pc/pc, (A19) 

where these equations are only dependent on the cooling function 
via the thermal state at which cooling vanishes, po/ po - Eliminating 
Pc pc gives 

u 2 s po - u s vopo - po = 0. (A20) 

The solution for the shock velocity u s = vo/2 — 
\J (vo/2) 2 + po/ ' po- Pc , Pc can be found by substitution. 

The conditions immediately to the right of the shock (v s , p s , 
p s ) can be found from the usual Rankine-Hugoniot relations, 

(vq - u s )po = (v s - u s )p s (A21) 
(vq - u s ) 2 po + po = (v s - u s ) 2 p s + p s {All) 

1/ x2 . 7 PO l f ,2 . 7 Ps / a 

-(vo-u s ) H -— = -(vs-Us) + -— .(A23) 

2 7-I po 2 7 - 1 p s 



